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proposed, to improve the performance. A numerical evaluation of the methods, when 
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EXECUTIVE SUMMARY 


Emitter localization comprises an important aspect of communications and is 
valuable in military and civilian applications. The importance of position information of 
a wireless emitter can be visualized if we consider the emergency 911-call situation when 
a mobile phone is used. Also, there are military applications in which it’s desirable to 
localize a non-cooperative electromagnetic source. 

The most widely used radio frequency (RF) position localization (PL) technique 
that appears to be superior to others, is the hyperbolic position localization, also known as 
the Time Difference Of Arrival (TDOA) technique. The TDOA utilizes cross-correlation 
techniques to estimate the delay of a propagating signal received at two different 
receivers. This delay estimation defines a hyperbola of constant range difference from 
the receivers, which are located at their foci. The intersection of a set of three or more 
hyperbolas provides the PL estimate of the source. 

This thesis is focused on the TDOA estimation, applied to GSM signals, 
generated with the use of the HP-ADS software, in noisy channels. Improvements in 
denoising, in conjunction with wavelet processing, are proposed for estimating the 
TDOA of signals received at two spatially separated sensors. We propose two wavelet 
based denoising methods: the Modified Approximate Maximum Likelihood (MAML) 
method and the Modified Fourth Order Moment (MFOM) method. 

Wavelet decomposition is applied to the output of both receivers. The wavelet 
coefficients are weighted, in each subband by appropriate factors, which are based on the 
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received signal statistics. The weighting aims to enhance the frequency bands where the 
signal is strong and to attenuate the bands when the noise is dominant. The TDOA 
estimate is obtained, by extracting the lag at which the maximum of the cross-correlation 
function of the denoised data occurs. 

Several situations are studied, primarily the one when unequal signal-to-noise 
ratio (SNR) conditions prevail, in the presence of Gaussian noise. The performance of 
the proposed denoising methods in a jamming environment is also addressed. A 
frequency modulated (FM) noise serves as jamming waveform in the simulations, where 
we allow the presence of jamming in one of the two channels. The Mean Squared Error 
(MSE) of the estimated delay is used as criteria of goodness (i.e., low MSE denotes a 
small error and thus a good localization). 

Simulation results indicate good performance of the denoising methods and 
improved estimates relative to the ones obtained using no denoising. In particular, the 
application of the MAML method allows a 70.4 to 91.37 percent improvement in the 
TDOA estimation, depending on the sampling interval. The respective improvement, 
applying the MFOM method, ranges from 62 to 79 percent. Jamming degrades the 
performance but still the extracted estimates are improved. A higher degree of 
correlation of the signal components reduces the probability of error in TDOA 
estimation. 

In addition, we note that in all cases studied the MAML method outperforms the 
MFOM method. The MAML method appears to be more robust in presence of Gaussian 
noise and jamming, while the MFOM method shows an undesirable “sensitivity” at low 



SNR levels (i.e., -3 dB or less) and when jamming is present. Finally the limits (lowest 
SNR values), in each channels are obtained above which denoising is beneficial. 
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I. INTRODUCTION 


A. BACKGROUND 

Wireless mobile radio systems have generated many new services in recent years. 
Cellular communication systems provide data and voice communication services with the 
convenience of mobility to their users. With the development of new Personal 
Communication Systems (PCS), it is inevitable to expect an increase in the type of 
services and in the availability of them. 

However, one important service, that current wireless radio systems are unable to 
offer, is efficient and enhanced 911 emergency services. In June 1996, the Federal 
Communications Commission (FCC) adopted specific regulations, requiring wireless 
service providers to implement E-911 services [1]. These rules seek to improve the 
reliability of wireless 911 services and to provide emergency service personnel with 
location information that enable them to locate and offer assistance to E-911 callers [2]. 

In September 1999, the FCC revised its rules to better enable carriers to use 
handset-based location technologies to meet the 1996 requirements [3]. In particular, 
FCC established separate accuracy requirements and deployment schedules for network- 
based and handset-based technologies. In brief the localization standards are: 

• Handset-based solutions: 50 meters for 67% of calls and 150 meters for 95% 
of calls 

• Network-based solutions: 100 meters for 67% of calls and 300 meters for 95% 
of calls, 

both starting October 1, 2001. 
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B. OBJECTIVE 

Several different position location (PL) techniques are candidates for a wireless 
PL system [4]. Radio frequency PL systems dominate the field since they appear to be 
more advantageous, with regard to low cost, ease of integration and potentially high 
accuracy. Furthermore, these techniques work with no need of modification in the 
existing cellular/PCS infrastructure. 

The most widely used RF PL technique that appears to be superior to others, is the 
hyperbolic position localization, also known as the Time Difference Of Arrival (TDOA) 
technique. TDOA utilizes cross-correlation techniques to estimate the delay of a 
propagating signal received at two receivers. This delay measurement defines a 
hyperbola of constant range difference from the receivers, which are located at their foci. 
The intersection of a set of three or more hyperbolas provides the PL estimate of the 
source. 

In this thesis, we will discuss the TDOA estimation and its application to Global 
System of Mobile (GSM) signals. We will investigate several methods to improve the 
TDOA estimate and the performance of these methods when unequal SNR conditions 
prevail. Furthermore, the sensitivity to jamming will be studied. 


C. THESIS OUTLINE 

The first chapter gives a general background and an introduction to the problem. 
Chapter II introduces the GSM system. Chapter IE presents the various position 
localization techniques that are used, with a more detailed discussion of the TDOA and 
the extraction of its estimate. Chapter IV examines the wavelet based denoising methods, 
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to improve TDOA estimation. Chapter V discusses the GSM test signals and the 
jamming waveform generation used in the simulations. Chapter VI presents the unequal 
SNR and the jamming simulation results. Finally, Chapter VII contains the conclusions 
as well as suggestions for future work on this topic. 
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II. GLOBAL SYSTEM FOR MOBILE (GSM) 


A. GSM HISTORY 

The Global System for Mobile communications is a digital cellular 
communications system. It was developed in order to create a common European mobile 
telephone standard but it has been rapidly accepted worldwide. 

In 1982 the Conference of European Posts and Telecommunications (CEPT) 
formed a study group to analyze and develop a Pan-European public land mobile system 
[5]. This system had to meet certain criteria such as 

• Support of International roaming 

• Low terminal and service cost 

• Good speech quality 

• Spectral efficiency 

• ISDN compatibility (Integrated Services Digital Network) 

• Support a variety of new services and facility 

In 1989, the European Telecommunication Standards Institute (ETSI) took the 
responsibility and the phase I GSM specifications were established in 1990 [6]. 
Commercial use and service was started in 1991 and by 1994 there was a 1.3 million, 
subscriber population worldwide, which grew to almost 380 million by December 2000. 

In 1995 the phase II GSM specifications where completed, which included the 
adaptation of North America’s PCS 1900, which mainly aimed to standardize the 
coverage of rural areas. 

Today more than 200 GSM networks (including PCS 1900) are operational in 110 
countries around the world. The GSM Association estimates that by the end of 2001 the 
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number of subscribers will be around 500 million, and by the end of 2004, this number 
will increase to 730 million worldwide [6]. 

B. GSM RADIO SYSTEM 

GSM utilizes two bands of 25 MHz, one for the forward link (base station to 
mobile) at 935-960 MHz and the 890-915 MHz for the reverse link (mobile to base 
stations). The combined use of Time Division Multiple Access (TDMA) and Frequency 
Division Multiple Access (FDMA) schemes, provide the ability of simultaneous access 
for multiple users. 

The available bandwidth of 25 MHz is divided among 124 channels of 200 KHz 
each. It should be noted that there are 125 channels. One is not used on purpose to 
prevent interference from other neighboring non-GSM systems. This channel is divided 
into a 100 KHz guard band at the upper and the lower end of the GSM spectrum. 
Different carrier frequencies are used for each of the 124 channels. Each individual 
channel is time division multiple accessed by users at different locations within a cell 
site. Time is segmented into intervals called frames. Each frame is further partitioned 
into assignable, non-overlapping, user time slots. The frame duration is 4.615 
milliseconds and is divided among eight timeslots, each of which is occupied by an 
individual user. Radio transmissions on both the forward and the reverse link are made at 
a channel data rate of 270.833 Kbps, using binary Gaussian Minimum Shift Keying 
(GMSK) modulation with a 3dB bandwidth-bit-period product (BT) of 0.3. 
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This data rate implies a bit duration of 3.7 pscc and an effective channel 
transmission rate per user, of 33.854 Kbps (270.833 Kbps/8 users). 

The total number of available channels within the 25 MHz band is 1000 (125 x 8). 
The GSM frame structure is presented in Figure 2.1 [7]. 
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Figure 2.1. GSM Frame Structure “From Ref. [7].” 

It is also worth mentioning that there are higher order frames, called multi-frames, 
within the GSM structure. The highest order frame is called a hyper-frame and has a 
time duration of 3 hours 28 minutes and 52.76 seconds. This long period of a hyper¬ 
frame is required to support encryption with high security. The extended and detailed 
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description of the higher order frames are out of the scope of this thesis and the reader 


can find the details in [7, 8]. 


C. GMSK MODULATION 

The Gaussian Minimum Shift Keying (GMSK) is a binary modulation scheme, 
which can be viewed as a derivative of the Minimum Shift Keying (MSK). In GMSK, 
applying a baseband Gaussian-shaped filtering to the Non-Retum-to-Zero (NRZ) 
modulated data prior to modulation reduces both the side lobe power level and the width 
of the main lobe. Figure 2.2 illustrates a simple block diagram of a GMSK modulator. 
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Figure 2.2. Simple GMSK Modulator 

The use of the pre-modulation Gaussian filtering smoothes the phase trajectory 
and hence stabilizes the instantaneous frequency variations over time, which leads to the 
reduction of the side lobe power spectral levels [8]. However, intersymbol interference 
(ISI) is introduced, because filtering converts the full response message signal into a 
partial response scheme where each transmitted symbol spans several bit periods T [9]. 
This is evident if we examine Figure 2.3, which depicts the impulse response of a 
Gaussian filter for a BT product of 0.3 and 0.5. 
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Figure 2.3. Gaussian Filter Impulse Response for BT = 0.3 and BT = 0.5 

Referring to Figure 2.3, we notice that a bit is spread over approximately 3 bit 
periods for BT = 0.3 and two bit periods for BT= 0.5. This is the reason that ISI 
phenomena arise with the use of the pre-modulation Gaussian filter. It is clear that as the 
product BT decreases, the ISI increases. 

A plot of the spurious radiated power in the adjacent channel to the desired 
channel power, with normalized frequency separation (f-f s )T (for GSM is 

3 2 

200-10 --= 1 . 5 ) on the abscissa and the 3dB bandwidth-bit-duration product 

270 . 833 - 10 3 

BT of the pre-modulation Gaussian filter, as parameter, is shown in Figure 2.4 [8]. 
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Figure 2.4. Power Spectral Density of GMSK (BT as parameter) “From Ref. [8]” 

We can see that while the GMSK spectrum becomes more compact (Figure 2.4) 
as the BT decreases, the degradation due to ISI increases (Figure 2.3). GMSK sacrifices 
the irreducible error rate, due to ISI, for very good spectral efficiency due to the constant 
envelope properties. 

For the GSM systems, GMSK with a BT of 0.3 at a channel data rate of 270.8 

Kbps has been adopted, since it appears to be an optimum compromise between bit error 
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rate and out-of-band interference. The value 0.3 represents the 3 dB bandwidth-bit- 
duration product or normalized pre-Gaussian bandwidth, which corresponds to a filter 
bandwidth of 81.25 KHz for a bit rate of 270.833 Kbps. 
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III. POSITION LOCATION TECHNIQUES 


A. TECHNIQUES 

It is inevitable that wireless communication carriers must have the ability to 
provide localization of every mobile device. This inevitability is validated by the Federal 
Communications Commission (FCC) requirement for the wireless emergency 911 
services [1, 2, 3, and 4], in the United States. 

In order to comply with the FCC ruling the carrier must provide positioning 
solutions that do not require the users to take any action in order to be located. There are 
three basic approaches to obtain position information that do not require modification of 
the mobile units, which are briefly described below. 

1. Angle of Arrival (AOA) 

The Angle of Arrival method utilizes multi-array antennas at cell site locations 
and tries to estimate the direction of arrival of the signal of interest. This implies that a 
single AOA measurement restricts the source location along the line in our estimated 
direction. In general multiple AOA estimates are needed to improve the unit’s location 
estimate, which is based on the point of intersection of the projected lines drawn out from 
the cell site at the angle from which the signal originated. 

Algorithms that exploit the phase differences between closely spaced antenna 
elements of an antenna-array provide the estimate of AOA. The spacing between the 
elements of the array is typically less than Vi wavelength of the received signal. This is 
required to avoid ambiguity in the estimating of the arrival angles. The resolution of 
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AOA estimates improves as the distance between the elements increases, however this 
improvement is at the expense of introducing ambiguities. 

Although angle of arrival has been well developed among military and 
government organizations, since it requires no special modifications to mobile devices, 
there are certain drawbacks. 

One important issue is that AOA is extremely sensitive to wide-angle reflections 
(multi-path) that occur when a device is operating in heavily shadowed channels, such as 
those encountered in urban environments. These reflections may provide power levels 
stronger than the Line-of-Sight (LOS) direction, thus an incorrect calculation of angle. 
However, the application of a technique known as Correlation Interferometry Direction 
Finding (CIDF) [10] provides efficient solution to the problem. This technique employs 
cross-correlation and uses a calibration database to determine AOA. The technique 
provides automatic compensation for the errors that the multi-path may introduce. The 
effectiveness of this method is evident from the fact that it has been successfully applied 
to shipboard High Frequency Direction Finding (HFDF) system, as discussed in 
Reference [11]. 

Another factor is the considerable cost of installation of antenna arrays, although 
the future use of “smart antennas” [12], which are compact arrays that shape the cell site 
transmitter and receiver signals into a beam that focuses on the mobile unit, is very 
promising. 

Finally the accuracy of the AOA estimate is reduced as the mobile unit moves 
away from a cell site or moves from cell to cell and the call is handed off from channel to 
channel. 
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2. Time Of Arrival (TOA) 

Time of arrival is another well-known technique used to determine the location of 
mobile units. The key point in this method is to determine the time it takes, for the 
signal, to travel from the source to the receiver on the forward or the reverse link. For 
this purpose, the base station transmits an inquiry to the mobile unit and measures the 
time in which the unit responds. The total measured time corresponds to the round trip 
signal delay plus the processing delay of the mobile unit. 

If the processing delay is known with sufficient accuracy then the round trip delay 
can be determined. Thus, using this time measurement, the distance between the unit and 
the receiver can be estimated. Taking the estimated distances from other two base 
stations, the position of the mobile is extracted by triangulation. 

It is obvious that timing errors in the absence of Line-of-Sight (LOS) are 
significant and very important in the estimation of the unit’s location. The impact of the 
multi-path reflections appears to be very strong, decreasing the accuracy of the estimated 
position. In practice it might be difficult for the base stations to estimate the response 
delay of the mobile unit with sufficient accuracy. 

3. Time Difference Of Arrival (TDOA) 

This method is well suited to estimate the location of wireless devices, since it is 
applicable with both brief transmissions, such as the reverse control channel, and with 
longer transmissions, such as the reverse voice channel [13]. This technique is based, in 
principal, on estimating the difference of arrival times of a signal received at multiple 
receivers. 
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The most common approach to obtain TDOA estimates, applies cross-correlation 
between the signals arriving at two spatially separated base stations. The lag that 
maximizes the cross-correlation functions represents the TDOA estimate. 

In addition, the cross-correlation is used to determine at which base station the 
signal arrived first. The combination of the estimates of two or more pairs of base 
stations yields hyperbolic curves along which the mobile may lie. The intersection of 
these curves provides the position estimate. 

We have to point out that it is essential, for the accuracy of this technique, that all 
of the base stations have a synchronized time reference. This means that a precise timing 
of a mobile’s signal transmission must be maintained at various cell sites. This is 
accomplished by using a common time base at the base stations. 

This technique offers many positive benefits as compared to other position 
localization techniques. Although both AOA and TDOA techniques rely mainly on a 
direct LOS between mobile and base stations, multi-path components affect TDOA less 
[13]. This is true in situations where the multi-path reflections affect all the stations 
involved in the localization process. The timing errors introduced may be cancelled or 
reduced in the time difference process [14]. Additionally TDOA does not require 
knowledge of the absolute time of transmission from the unit, which in general is more 
difficult to determine. Lower implementation cost is also an advantage since TDOA 
techniques can be used with the existing equipment and simple, inexpensive conventional 
antennas (i.e., whip antennas). 


16 



B. TDOA ESTIMATION AND POSITION LOCALIZATION 


The position localization process consists of two distinct parts. One is the 
computation of accurate TDOA estimates. This process in general involves denoising 
techniques that aim to discard the undesirable presence of noise in the signal of interest. 
The most popular denoising techniques are based on wavelet decomposition and will be 
presented in detail in the next chapter. 

The second part of the process uses the TDOA estimates to extract the position of 
the unit. We will briefly discuss both processing stages and point out the principle idea 
behind each one. 

1. Determination Of TDOA Estimates 

The most popular method to obtain TDOA estimates is the generalized cross¬ 
correlation method. The transmitted signal s(t), when it is received by the base stations 
will be delayed and corrupted by noise and can be represented as follows 

x (t) = a s (t - d x ) + n x (t) (3.l) 

y{t) = fi s(t-d y )+n y {t) (3.2) 

where d x , d y are the time delay from the transmitter unit to receiver x and y respectively 
and, a, /3 are coefficients such that 0 < a, J3 <1. 

Bearing in mind that the receiving base stations are synchronized, so that a 
common time base is established, the cross-correlation function of the received signal is 
computed as [13] 
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K ( T ) = b \ x (*) y( t + T ) dt ( 3 - 3 ) 

1 0 

Once R (t) is computed, the lag that maximizes (3.3) is considered to be the estimate 
for the TDOA value. 

2. Determination Of Position Estimate 

When the previous process is completed, the TDOA estimates are available. Each 
estimate localizes the unit on a hyperboloid with a constant range difference between the 
two base stations. This range difference is given by 

(1 ',-yf+ (z,-z) ! -j(x J -xf+ (y,-T+ (z,-3 2 (3.4) 

where (X i , Y i , Z,.) and (X j , Y j , Z y ) represent the coordinates of the base stations i and 

j and (x, y, z) are the possible mobile unit coordinates [13]. Of course, the mobile unit 
only occupies one point on the hyperbolic curve. 

Using three pairs of base stations and substituting the TDOA estimates in 
equation (3.4), we can estimate the location of the mobile unit. Several methods are 
available for solving the hyperboloid equations with great accuracy and a detailed 
presentation can be found in References [15, 16, 17]. 
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IV. DENOSING METHODS TO IMPROVE TDOA ESTIMATE 


In this chapter we will present the details of the denoising methods, used to 
improve the accuracy of the TDOA estimate. This is a crucial issue, that directly 
influences the position localization estimate of a mobile unit, as pointed out in section 

m.B.2. 

Generally speaking, the denoising process consists of three phases: (a) wavelet 
transform of the received signal, (b) thresholding of the decomposition coefficients and 
(c) inverse wavelet transform of the modified coefficients. 

A signal emanating from a remote source, embedded in noise and sampled at two 
spatially separated sensors can be modeled as 

x (k) = s (k) + n x (k) (4.1) 

y{k) = s{k + d) + n y {k) k = 0,1,2,..., A -1 (4.2) 

where s(k) is the unknown signal, which is assumed to be independent of the noise and a 
and /? coefficients are set to one. The components n x (k) and n y (k) are additive 
uncorrelated Gaussian noise sequences that control the signal-to-noise ratio (SNR). The 
quantity cl is the time difference of arrival and N is the number of samples collected at 
each receiver. 

The methods proposed to increase the accuracy of the TDOA estimation are 
generically depicted in Figure 4.1. The noise at each receiver is effectively reduced, by 
using wavelet based denoising on both channels. 
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Wavelet Denoising 


Correlator 


Figure 4.1. Block diagram for TDOA estimation using wavelet denoising 


A. MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD (MAML) 
METHOD 

The Approximate Maximum Likelihood (AML) method was proposed by Y.T. 
Chan, H.C. So, and P.C. Ching for estimation of the TDOA of signals [18]. It aims to 
attenuate the noise in the subbands where the noise is strong and to enhance the signal at 
the subbands where it is strong, by appropriate weighting of the decomposition 
coefficients. The prefilters H x (f) and H y (f), shown in Figure 4.2, are chosen to satisfy 
[18] 




Gjf) 

, G„(/) ! Gj/T 

G.Jf) 


(4.3) 
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where G ss (f), Gn x n x (f), and Gn y n y (f) denote the auto-power spectra of s(k), n x (k), and n y (k), 
respectively. This choice of the pre-filters is known as maximum likelihood (ML) 
weighting. 



Figure 4.2. Generalized cross-correlator 


The AML filtering is performed by weighting the noisy signal in each 
decomposition level. The weights are determined level by level, so that a large delay 
variance is avoided. The derivation of these weights assumes that the signal component 
in each subband can be approximated by a flat spectrum. This assumption leads to the 
idea that the optimal weighting function and the auto-power spectrum of the signal, are 
modeled by a piecewise constant functions. This makes sense, particularly when the 
number of the decomposition levels is sufficiently large. 

Using the above, the weights for the detail coefficients are given by 


w 


d, 


sdj 


+ <7 


Sd: 


+ <7; 


(4.4) 
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where, i=l,2,...J, and /is the number of the wavelet decomposition levels, [18]. 

Under the assumption that the wavelet filter’s gain is two and the filter transfer 
function is ideal, the noise power will remain the same after passing through each filter. 
The estimates for the noise and the signal powers are computed as follows: 


arg max 


N -1 


A k=0 


< 7„ 2 


1 N -1 

K =^ 2>®) 2 -*, 2 

iy k =0 


1 


T:—l 


=jYX d -A k jf ~ d l 


k =0 


(4.5) 

(4.6) 

(4.7) 

(4.8) 


The quantity 6 2 s denotes the estimate of the total input signal power, d 2 and d] 

represent the estimated noise powers, and d 2 d is the estimated signal power in the 

subband The quantities N and T, represent the number of elements of the signal 
sequence and the number of elements in each subband, respectively. A modified weight, 
compared to the one proposed in [19], is used for the last approximate coefficient of the 
wavelet decomposition, as follows: 
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w 


(4.9) 


G 2 

_ _ SOj _ 

a J — *2 . ~ 2 TTl , ~ 2 \ 

G G -h G I G + G I 

n l n 2 saj V n x n 2 / 


where 

j 

*2 r\ J *2 X 1 r\ J—i *2 

&SU, = 2 2 ^ 

1=1 


(4.10) 


The d 2 sa is the estimated signal power in the approximate coefficient of the wavelet 
decomposition at level J. The modified coefficients are obtained as: 

d, = w d ■d i i = l,2 ,...J (4.11) 

‘mod “/ * 


After weighting each decomposition coefficient, the inverse wavelet transform is 
used to reconstruct the denoised signal, from the modified coefficients. The adaptation of 
the AML method and its application in both channels is referred to as the Modified 
Approximate Maximum L ikelihood (MAML) method. 

One important issue is to determine the number of the appropriate decomposition 
levels, in order to avoid considerable variances in the power estimates (especially when 
the number of elements in the decomposition coefficient is small). A simple assumption 
is made to avoid large variances. If the adjacent subband signal power is small then we 
interrupt the decomposition. By trial and error, for the GSM signal that we are using, we 
found that if the ratio of the power of the approximate coefficient, to the power of the 
detail coefficient, in a given scale (level), is smaller that 0.5, we can terminate the 
decomposition procedure without a measurable loss. 
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B. MODIFIED FOURTH ORDER MOMENT (MFOM) METHOD 


1. Background 

White Gaussian noise has samples that are statistically independent. If z is a 
random vector of N consecutive samples from a real valued, zero mean, white Gaussian 
process, its probability function is given by [20]: 


N-l i 

/(z 0 ,z 1 ,...,z iV _ 1 )= JJ-==-c 2<T °" (4.12) 

n= 0 ^)27TCT ( j 

The moments of a stationary process are given by [20] 

= E{z[n§=m z (4.13) 

M ! 2) [/, ] = E{z[n]z[n + /,]}= R z [/, ] (4.14) 

M | 3) [/,, l 2 ] = E{z[n]z[n +1 , ]z[n + l 2 ]} (4.15) 


M * 4) [l 3 , l 2 , l 3 ] = £{z[tt]z[« + l v ]z[n +1 2 ]z[n + l 3 ]} 

= £{z[n]z[n + /, ]}e{z[« + 1 2 ]z[n +l 3 ]} 

+ £{z[n]z|n + 1 2 ]}e{z[« + /j ]z[n +1 3 ]} 

+ £{z[tt]z[« + 1 3 ]}e{z[« + h ]z[« + U_ ]} (4.16) 

and 

m‘ 4) [0,0,0] = 3-(J 4 (4.17) 


The fourth order moment of the received signal, i.e. x(k)=s(k)+nj(k), is computed 
as follows: 
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(4.18) 


m[ 4) [oao] = 4(^)+^)) 4 } 

= E{s 4 [k]}+ 4E{s 3 [k]n l [k]}+ 6 e{s 2 [k]nf [/.]} 

+ 4£'{s[/c]rt | i [A:]}+ E{n*[k]} 

^ 3 < 

Going a little bit further, one can compute the mean and the standard derivation of 
the fourth order moment. If we define p = M[ 4 \ then the mean of p is: 

7li £ fe, 4 }=3<7* (4.19) 

N i=0 

In addition, the second moment and the variance are: 

8 

E{p 2 }= — ~(9N + 96) (4.20) 


a 


2 

P 



(4.21) 


2. Method 

In the proposed method we accounted for the variability of the estimate by 
choosing a threshold of 3. Id 4 , equation (4.19). 


It is evident that the fourth order moment of the detail decomposition coefficient, 
which contains signal, will be greater than 3d 2 , , where <j 2 nd is computed in equations 

(4.6) and (4.7) and denotes the decomposition level. This fact is used to eliminate the 
wavelet coefficients in which the noise is dominant and keep the ones where the signal is 
strong. The estimates of the signal and noise powers are computed according to 
equations (4.5), (4.6), and (4.7). 
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The noise power, after passing through the first high pass filter is 


<T„ 2 . 


<7; 


j = U 


(4.22) 


where G is the gain of the high and low pass filters, at decomposition level 1. 

Assuming that the filter gain is 2, then the noise power of any detail 
decomposition coefficient d, is 


G l dL. I 

■ <y„ 


2 ' 


j = 1,2 


(4.23) 


We modify each wavelet decomposition coefficient as follows: 

4 ; ifEfc 4 } a 3.1-< 

0 ; otherwise 



(4.24) 


We introduce a modification in the above presented method. Prior to the 
application of equation (4.24) an additional threshold is used to pre-filter the wavelet 
decomposition coefficients in time. Based on the previous discussion, the power 
contained in each element of the decomposition coefficient, which represent signal, will 

be greater than ^3.1(7^ . Thus, the empirical threshold corresponding to the square root 
of 3.1(7^ is used. This threshold is compared with the squared value of each element of 

the decomposition coefficients. If the elements have a squared value less than ^3.1 61 

they are set to zero, otherwise they are left unchanged. After modifying the coefficients, 
the fourth order method described previously is applied. Finally, the resultant modified 
coefficients are combined to reconstruct the denoised signal using the inverse wavelet 
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transform. We have to mention that the procedure described in Section IV.A (i.e. 


minimize the number of decomposition levels) is also applied. 


C. COMBINED DENOISING METHODS 

In our effort to eliminate the noise in the GSM noisy signal, we try combining 
denoising schemes, using the previously presented methods. The combined schemes are 
depicted in the block diagrams of Figure 4.3. Even though the resulting denoised signal 
is satisfactory, the combined schemes do not outperform the individually applied 
denoising methods. 

The results and the performance of the combined schemes, as compared to the 
individually applied over, are presented in Appendix [A]. 


Noisy 

Signal 


(a) 



Denoised 

Signal 


Denoising 


Noisy 

Signal 


(b) 



Denoised 

Signal 


Denoising 

Figure 4.3. Combined denoising schemes 
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V. TEST GSM SIGNALS AND JAMMING SIGNAL DESCRIPTIONS 


A. GSM SIGNAL 

The Hewlett-Packard Advanced Design System (HP-ADS) was used to generate 
the test GSM signals, for our simulation purposes. 

The HP-ADS is a very powerful and reliable tool for engineers [21], since it 
offers the ability to the user to customize and evaluate many different design aspects of 
communications, digital signal processing, and a variety of other system parameters. The 
HP-ADS software provides an excellent method to obtain signal details that are 
essentially identical to the ones encountered in an actual cellular communication system. 
Detailed operating instructions of the HP-ADS program and a more specific view of each 
of its components can be found in Reference [22]. 

Using the communications package of the HP-ADS, we generate a GSM signal, 
sampled at three different sampling rates [22]. The specifications of the GSM signal 
discussed in Chapter II are used as inputs. The specified symbol period of 3.7 
microseconds is used. The filter bandwidth at the receiver is 1.2 MHz, which implies that 
the minimum, non-violating the Nyquist criterion, sampling frequency is 2.4 MHz. The 
parameters of the test signals are summarized in Table 5.1. 
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Samples/Symbol 

Sampling Interval 

(nsec) 

Sampling 

Frequency (MHz) 

Symbols per 

600 samples 

40 

92.59 

10.8 

15 

20 

185.185 

5.4 

30 

10 

370.370 

2.7 

60 


Table 5.1. Test GSM signal parameters 


The I and Q channel outputs of each test signal, are shown in Figures 5.1, 5.2, and 
5.3. It is evident that as the sampling interval increases, so does the variation in the GSM 
signals. 

The reason for this is obvious, if we consider that for a constant number of 
samples (600), we are using in all simulations, a lower sampling frequency captures more 
information (symbols) than a higher one. 

(a) I Channel (b) O Channel 




Number of samples Number of samples 


Figure 5.1. I and Q channels of the HP-ADS data with T s =92.59 (nsec) 
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(a) I Channel 


(b) Q Channel 




Number of samples Number of samples 


Figure 5.2. I and Q channels of the HP-ADS data with T s =185.185 (nsec) 


(a) I Channel (b) Q Channel 




Number of samples Number of samples 

Figure 5.3. I and Q channels of the HP-ADS data with T s =370.370 (nsec) 

At this point, we have to mention the auto-correlation characteristics of this set of 
samples. The respective auto-correlations of the test signals are shown in Figure 5.4. 
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(a) HP-ADS data with T g =92.592 (nsec) (b) HP-ADS data with T s =185.185 (nsec) 




(c) HP-ADS data with T s =370.370 (nsec) 



Figure 5.4. Auto-correlation of the HP-ADS GSM test signals 
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Note that a better correlation is observed with increase in the sampling interval. 
The main lobe is sharper, while the side lobes are reduced. This is a key point in the 
upcoming discussion regarding the results of our simulations. 

B. JAMMING SIGNAL 

In order to demonstrate the performance of the proposed denoising methods in a 
jamming environment, a noise jamming waveform is used. 

The purpose of the noise jamming is to interfere with the enemy’s electronic 
equipment so that the actual signal is distorted or over shadowed. The advantage of using 
noise jamming is that very little detail of the enemy’s equipment needs to be known. 

Based on information theory, the noise with a uniform spectral density, over the 
band of interest, provides one of the best noise-jamming waveforms. That is the white 
Gaussian noise, which has the maximum entropy, compared to other random waveforms, 
for a specified average power [23]. This characteristic makes it difficult to distinguish 
between receiver noise, in the victims system, and externally injected jamming noise. 

The simplest method to generate a noise-jamming signal is high-power 
amplification of band limited Gaussian noise. This method is known as direct noise 
amplification (DINA) [23]. 

However, most of the high-power microwave amplifiers are limited in the peak 
power they can handle. For this reason, alternative and much more flexible methods have 
been established to generate high power noise. One of these utilizes frequency or phase 
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modulation of the jammer’s signal, with a random waveform. The advantage of this 
method is that maximum output power can be obtained from the jammer transmitter. 

Based on the above, we are using a jamming signal referred to as “erfed” FM 
noise [23]. A block diagram, shown in Figure 5.5, illustrates the principles of generating 
this type of noise to be used as jamming waveform. 


Gaussian 

noise 

source 



Erfed 

signal 


Figure 5.5. Block diagram of Erfed FM noise generation 

In Figure 5.5, a base band Gaussian frequency control waveform (noise) is 
applied to frequency modulate a power oscillator, such as a backward wave oscillator 
(BWO) or carcinotron microwave source. This allows the buildup of a Gaussian pseudo 
noise jamming waveform in the victim’s receiver. In order to produce the desired 
uniform RF power spectrum at the output of the jammer, it is necessary to alter the 
amplitude distribution of the modulating noise before it is used to modulate the carrier 
[24]. It has been shown that the necessary transformation can be achieved with the use of 
a nonlinear network with transfer function the error function (erf()) [24]. That is why, the 
output waveform from the VCO, is called “Erfed” noise. 

A MATFAB code is used to generate the desired jamming waveform. Two 

scenarios are simulated. One generates a jamming signal with a bandwidth equal to the 

bandwidth of the receiver (1.2 MHz) [22], while in the second; the jamming waveform 
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interferes only with the bandwidth of a physical channel, which is 200 KHz. The base 
band Gaussian noise source used in these scenarios has a bandwidth of 55 KHz, for the 
first scenario, and a 20 KHz, for the later. 

Figures 5.6 (a) and 5.7 (a) show the Gaussian noise waveform for the two 
scenarios, while Figures 5.6 (b) and 5.7 (b) depicts their respective Power Spectral 
Densities (PSD). Figures 5.6 (c) and 5.7 (c) display the generated noise spectra for the 
two test cases, while Figures 5.6 (d) and 5.7 (d) show the time behavior of the jamming 
signal. 
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(a) Gaussian Freq. control waveform 


(b) Spectrum of the control waveform 



Figure 5.6. 


Erfed Jamming waveform with bandwidth 1.2 MHz 
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(a) Gaussian Freq. control waveform 


n 


(b) Spectrum of the cpntrol waveform 



DO 

■o 



Hz 


x 10 


(d) Jamming waveform 



Figure 5.7. Erfed Jamming waveform with bandwidth 200 KHz 


The detailed MATLAB code for the generation of the jamming waveforms is 
shown in Appendix [C]. 
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VI. SIMULATION RESULTS 


A. UNEQUAL SNR CONDITIONS 

In this section we apply the methods described in Chapter IV, to our test HP-ADS 
generated signals embedded in zero mean, white Gaussian noise. 

Our goal is to investigate how well we can extract useful TDOA data. The Mean 
Squared Error (MSE) is used as the measure of improvement. We define the MSE as the 
square of the difference between he true TDOA value and the estimated TDOA value. 
We compare the MSE obtained from the proposed methods (improved value), to MSE for 
non-denoised TDOA estimation (original value). 

Each of the proposed denoising methods is evaluated for each of the three 
sampling frequencies. For the convenience of the reader we label the test signals as 
follows: 

(1) HP-ADS signal with sampling interval of 92.592 nsec. 

(2) HP-ADS signal with sampling interval of 185.185 nsec. 

(3) HP-ADS signal with sampling interval of 370.370 nsec. 

All simulations were performed using MATLAB [25]. The codes used are 
provided in Appendix [C]. For lucidity in the following discussion, we define the 
percentage improvement and the average percentage improvement. The percentage 
improvement is given by 

MSE value without denoising - MSE value with denoising 

Percentage improvement = - x 100 

MSE value without denoising 
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The average percentage improvement is simply the mean value of the percentage 
improvement for all values of total SNR. True TDOA is randomly chosen between -150 
to +150 samples, thus the MSE value represents the sample difference squared. To 
quantify the MSE further, consider for example a value of 4. This value corresponds to a 
sample difference of 2 and hence a time delay of 2x(Sampling Interval (nsec)). The 
transforms use a Daubechies wavelet of order 8 (i.e., db4 [26, 27]), and the optimum 
level of decomposition was determined, by trial and error, to be 5. Unequal levels of 
noise are used in each channel, and the lowest SNR value (limit) is obtained, that still 
permits reasonable localization. Reasonable localization provides a MSE smaller than 
that obtained by using raw data. A fixed value of SNR ranging from 3 dB down to -9 
dB, is kept in one channel while the SNR values in the other channel are varied [28]. The 
signal power is unity and different SNRs, are obtained by scaling the random noise 
sequences. The MSE versus the total SNR in both channels is plotted. The total SNR 
represents the sum of the SNR values of the two channels. The results are shown in 
Figures 6.1, 6.2 and 6.3. Figures 6.1(a), 6.2(a), and 6.3(a) depict the region of the SNRs 
over which we have results that permit reasonable localization, while Figures 6.1(b), 
6.2(b) and 6.3(b) show the limiting cases (lowest SNR in each channel). 

One obvious trend is that the MSE improves, i.e., decreases, as the sampling 
interval increases. This observation is also evident in Tables 6.1 and 6.2, which contain 
the numerical values of the computed MSE for each test signal as obtained by using the 
MAML method and the Modified Fourth Order Moment (MFOM) method, respectively. 
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MSE MSE 


In both tables, we can clearly see a decrease in error values as we move across 
each row from left to right. Note that in Table 6.1 and for the MAML case for data with 
sampling interval 185.185 (nsec), a small inconsistency exists. 


(a) 



(b) 



-15 -10 -5 0 5 10 

Total SNR (dB) 


Figure 6.1. MSE versus total SNR for HP-ADS data with T s =92.59 (nsec) 
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(a) 



(b) 



Figure 6.2. MSE versus total SNR for HP-ADS data with T s =185.185 (nsec) 
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(a) 



(b) 



Figure 6.3. MSE versus total SNR for HP-ADS data with T s =370.370 (nsec) 


43 


















Mean Squared Error (MSE) 


Total 

>NR in 

Sampling Interval 

92.59 (nsec) 

Sampling Interval 

185.185 (nsec) 

Sampling Interval 

370.370 (nsec) 

(dB) 

Non Dcnoised 

MAML 

Non Dcnoised 

MAML 

Non Dcnoised 

MAML 

10 

8.20 

3.00 

1.83 

0.71 

0.67 

0.17 

8 

10.92 

3.69 

2.65 

0.71 

0.79 

0.17 

6 

15.92 

5.10 

3.50 

0.90 

0.88 

0.18 

4 

19.03 

8.04 

4.05 

1.10 

0.98 

0.22 

2 

24.78 

10.33 

5.73 

1.79 

1.22 

0.45 

0 

28.07 

10.95 

6.50 

1.78 

1.52 

0.46 

-2 

31.78 

12.74 

8.36 

1.92 

1.83 

0.51 

-4 

39.56 

15.51 

13.41 

2.31 

2.60 

0.61 

-6 

46.79 

20.04 

16.15 

3.16 

3.62 

0.80 

-8 

68.81 

25.50 

18.98 

3.76 

4.69 

1.00 

-10 

87.38 

33.07 

22.38 

4.64 

5.72 

1.27 

-12 

118.18 

44.19 

30.52 

6.13 

6.93 

1.72 

-14 



35.80 

7.33 

10.19 

2.63 

-16 



46.22 

11.11 

174 

4.38 

-18 



297 

17.87 

492 

12.25 


Table 6.1. MSE values at different values of total SNR for MAML method 




Mean Squared Error (MSE) 


Total 

>NR in 

Sampling Interval 

92.59 (nsec) 

Sampling Interval 

185.185 (nsec) 

Sampling Interval 

370.370 (nsec) 

(dB) 

Non Dcnoised 

MFOM 

Non Dcnoised 

MFOM 

Non Dcnoised 

MFOM 

10 

8.20 

3.52 

1.83 

0.65 

0.67 

0.19 

8 

10.92 

4.07 

2.65 

0.81 

0.79 

0.22 

6 

15.92 

5.68 

3.50 

1.07 

0.88 

0.27 

4 

19.03 

7.05 

4.05 

1.51 

0.98 

0.36 

2 

24.78 

8.59 

5.73 

2.18 

1.22 

0.63 

0 

28.07 

14.03 

6.50 

2.62 

1.52 

0.73 

-2 

31.78 

18.95 

8.36 

3.12 

1.83 

0.76 

-4 

39.56 

26.79 

13.41 

4.10 

2.60 

1.10 

-6 

46.79 

33.41 

16.15 

5.76 

3.62 

1.60 

-8 

68.81 

43.08 

18.98 

7.92 

4.69 

2.06 

-10 

87.38 

59.87 

22.38 

10.66 

5.72 

2.65 

-12 

118.18 

82.68 

30.52 

15.01 

6.93 

3.43 

-14 



35.80 

19.62 

10.19 

4.55 

-16 



46.22 

28.91 

174 

7.75 

-18 



297 

46.65 

492 

25.12 


Table 6.2. MSE values at different values of total SNR for the Modified Fourth Order 

Moment method 




This is the value of 1.79 observed at 2 dB total SNR, which is greater than the 1.78 
observed at 0 dB of total SNR. This may be due to the random number generator we use 
to generate the noise sequences, or to the number of the realizations (100 realizations) for 
the simulations purposes. An increase in the number of the simulation repetitions, may 
remove this inconsistency. 

The limits obtained for each method and for each sampling rate are given in Table 6.3. 


Sampling 

Interval (nsec) 

Lowest SNR allowed in each 

channel (dB) 

Maximum difference in SNR 

between the two channels (dB) 

MAML 

Modified 4 th 

Order Moment 

MAML 

Modified 4 th 

Order Moment 

92.59 

-8 

-6 

14 

12 

185.185 

-9 

-8 

18 

17 

370.370 

-9 

-9 

20 

18 


Table 6.3. Limitations in SNR values 


The limiting cases represent the threshold below that the denoising methods do 
not perform well. The MSE in these cases, appear to have larger values than the non- 
denoised MSE values, as computed by the direct cross-correlation of the raw data. 

An important fact is that in addition to the limits of the lowest SNR in each 
channel, there is also a restriction regarding the combination of the SNRs in the channels, 
that is the total SNR. For example, consider a value of -3 dB for total SNR. This can be 
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achieved by a large number of combined values of SNRs in each channel. But not all of 
the combinations are allowed; only the ones that have a dB separation smaller than the 
maximum difference. This maximum difference is found to be as shown in the last two 
columns of Table 6.3. Again we see that this limit increases as the sampling interval 
increases. The MAML method appears to permit a larger number of SNR combinations 
to achieve a specific total SNR value compared to the MFOM method. 

Interpreting Figure 6.1(a), we see that using the HP-ADS data with a sampling 
interval of 92.59 (nsec), the average percentage improvement over all total SNR values 
between -12 and 10 dB is about 61.53% for the MAML method, and 46.71% for the 
MFOM method. 

The HP-ADS data with a sampling interval 185.185 (nsec), see Figure 6.2(a), 
permits an average percentage improvement, over all total SNRs between -18 and 10 dB, 
of 76% when MAML denoising is applied, and 60.5% when the Modified Fourth Order 
Moment denoising is used. 

The average percentage improvement for the HP-ADS data with a sampling 
interval of 370.37 (nsec), see Figure 6.3(a), is even better. In this case we obtain 
improvements of 78% and 63.3% for the MAML and the Modified Fourth Order 
Moment, respectively, over total SNR values between -18 and 10 dB. The MSE for all 
methods and for all SNRs above a total of lOdB is essentially zero, provided each 
channel has an SNR above the lowest SNR value shown in Table 6.3. 


47 



The improvements we discussed above can also be related to the plots in Figure 
5.4, which presents the auto-correlation function of the test signals. There, we observed 
improved correlation with decreasing sampling frequency. 

Hence we can state that a higher degree of correlation of the signal components 
reduces the probability of error in the TDOA estimation. 

B. PERFORMANCE IN A JAMMING ENVIRONMENT 

In Chapter IV, we discussed in detail two denoising methods to improve TDOA 
estimates. In this section, we will use these methods to test the performance for the HP- 
ADS generated signals embedded in Gaussian white noise subject to a jamming signal, as 
described in Chapter V. 

Again, the goal is to investigate how well useful TDOA information can be 
extracted. The mean squared error (MSE) is the measure of improvement, i.e., the lower 
MSE the better the TDOA estimate. 

We define the percentage degradation of the MSE, and the average percentage 
degradation, in an analogous way as for the percentage improvement. The corresponding 
formula is 

„ ... (MSE value without Jamming)-( MSE value when Jamming present) _ 

Percentage degradation =--—-—- xlOO 

(MSE value without Jamming) 

Each denoising method is tried in two jamming scenarios, as described in Chapter 
V, for each of the three HP-ADS data sets. We will only present the case in which the 
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full bandwidth of the receiver is jammed, since the other case follows a similar trend. 
The results and the related plots of the latter case can be found in Appendix [B]. 

The labeling of the test signals is kept as in the previous section. For the jamming 
simulations, each channel uses equal SNR levels. The values of SNRs, used in both 
channels, range from 20 dB to -6 dB. Four values of jamming-to-signal (J/S) ratios are 
used, 10 dB, 5 dB, 0 dB, and -10 dB. These values are obtained by proper scaling of the 
jamming waveform. 

The MSE versus SNR in each channel is plotted, with the jamming-to-signal ratio 
as a constant parameter. The results are obtained using 100 realizations for each SNR, 
and are shown in Figures 6.5, 6.6, and 6.7. 

Figure 6.4 shows the performance of the denoising methods when no jamming is 
present. This provides a reference for how the denoising performs in a jamming free 
environment. Notice that the MSE improves (decreases), as the sampling frequency 
decreases and that the MAML method outperforms the Modified Fourth Order Moment 
method in all cases. 

The improvements can also be seen in Table 6.4, which contains the average 
percentage improvement computed for each method and each sampling interval. 
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(a) HP-ADS data with Ts=92.592 (nsec) 



(c) HP-ADS data with Ts=370.370 (nsec) 



Figure 6.4. MSE versus SNR for 


(b) HP-ADS data with Ts=185.185 (nsec) 



HP-ADS data, when jamming is not present 
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Average Percentage Improvement % 


Sampling Interval 

Sampling Interval 

Sampling Interval 


92.59 (nsec) 

185.185 (nsec) 

370.370 (nsec) 

MAML method 

70.4 

88.63 

91.37 

MFOM method 

62 

77.14 

79 


Table 6.4. Average percentage improvement in MSE, when jamming is not present 

The improvement in MSE is related to the sampling frequency and is related to 
the degree of correlation of the signal components. Improved correlation reduces the 
probability of error in the TDOA estimation. 

Continuing the discussion, lets take a look at Figures 6.5, 6.6, and 6.7. Each 
figure represents the results of denoising when jamming is present in one channel, for the 
HP-ADS data with sampling interval 370.370 (nsec), 185.185 (nsec), and 92.59 (nsec), 
respectively. 

We see that the MSE improves (decreases) as the sampling interval increases, or 
when the jamming-to-signal ratio decreases. It is obvious that higher levels of jamming 
power present in our signals degrade the performance of the denoising methods. The 
quality of MSE estimates degenerates as the power of the jammer increases. Also the 
SNR range, in which reasonable estimates are obtained, decreases as the J/S ratio 
increases. 
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(a) J/S ratio 10 dB (b) J/S ratio 5 dB 




(c) J/S ratio 0 dB 




-5 0 5 10 15 20 

SNR 


Figure 6.5. MSE versus SNR for HP-ADS data with T s =370.37 (nsec), at different J/S 

ratios 
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MSE MSE 


(a) J/S ratio 10 dB (b) J/Sratio5dB 




(c) J/S ratio 0 dB 



(d) J/S ratio-10 dB 



Figure 6.6. MSE versus SNR for HP-ADS data with T s =185.185 (nsec), at different 

J/S ratios 
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MSE MSE 


(a) J/S ratio 10 dB (b) J/S ratio 5 dB 




(c) J/S ratio 0 dB 



(d) J/S ratio-10 dB 



Figure 6.7. MSE versus SNR for HP-ADS data with T s =92.592 (nsec), at different J/S 

ratios 
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Specifically, using the HP-ADS data, with a sampling interval of 370.370 (nsec) 
(see Figure 6.5(a)) provides reasonable estimates above -3 dB for a lOdB jamming-to- 
signal ratio. Figures 6.5(b-d) allow reasonable MSE values for SNR’s greater than -5 
dB. Table 6.5 contains the average percentage degradation of the computed MSE for 
different values of J/S ratios, as compared to the MSE obtained when no jamming is 
present. 


Jamming-to-signal ratio 

(dB) 

Average Percentage Degradation % 

MAML method 

Modified 4 th Order Moment 

method 

10 

73.98 

88.71 

5 

70.05 

78.91 

0 

60.05 

65.24 

-10 

28.30 

33.97 


Table 6.5. Average percentage degradation in MSE at various J/S ratios (HP-ADS 
data with sampling interval 370.370 (nsec)) 


We also see that using the MAML method at a lOdB J/S ratio, the MSE values are 
on average, 73.98% worse than those obtained when no jamming is present. The MFOM 
method gives MSE values that are 88.71% worse for the same J/S ratio. It is evident that 
the MAML method outperforms the Modified Fourth Order Moment method, since it 
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gives MSE values which on average are better (relative to the no jamming situations) 
than those obtained using the MFOM method. 

For the HP-ADS data with a sampling interval of 185.185 (nsec) the results are 
shown on Figure 6.6(a-d) and the computed average percentage degradation is given in 
Table 6.6. 


Jamming-to-signal ratio 

(dB) 

Average Percentage Degradation % 

MAML method 

Modified 4 th Order Moment 

method 

10 

78.50 

93.68 

5 

76.74 

90.69 

0 

68.30 

81.60 

-10 

48.68 

54.36 


Table 6.6. Average percentage degradation in MSE at various J/S ratios (HP-ADS 
data with a sampling interval of 185.185 (nsec)) 


In this case, we observe that the presence of jamming limits the SNR value in the 
channels to be larger than -3 dB and in addition, we see an increased average percentage 
degradation, relative to the values of Table 6.5. 

The average percentage degradation at KMB J/S ratio is 78.5% using the MAMF 
method, and 93.68% when the MFOM method is applied. Still, the MAMF method gives 


56 




better results than the MFOM method, but the performance is reduced as the sampling 


frequency is increased. 

The HP-ADS data with a sampling interval 92.59 (nsec), see Figure 6.7(a-d), 
provides the worst MSE values. In this case, the presence of jamming severely affects 
the extraction of reasonable estimates. The SNR range is required to be larger than 6 dB 
when the J/S ratio is 10 dB, see Figure 6.7(a), while we require a 3 dB SNR in each 
channel when J/S ratio is 5 dB or less, see Figure 6.7(b). Table 6.7 shows the percentage 
degradation of MSE, based on the minimum SNR in each channel. 


Jamming-to-signal ratio 

(dB) 

Average Percentage Degradation % 

MAML method 

Modified 4 th Order Moment 

method 

10 

72.28 

98.78 

5 

69.68 

97.58 

0 

65.41 

94.75 

-10 

55.54 

80.17 


Table 6.7. Average percentage degradation in MSE at various J/S ratios (HP-ADS 

data with sampling interval 92.59 (nsec)) 


Considering the results presented in Tables 6.5, 6.6, and 6.7, we can say that the 
increase in sampling frequency reduces the quality of the TDOA estimation (compare the 
average degradation performance, under the same J/S ratio, of different sampling 


57 




intervals). In addition we observe that the MAML method out performs the MFOM 
method in all cases (in each Table MAML gives better results than the MFOM method 
does). The MAML method also appears to be more robust in a jamming environment, 
since the degradation of the computed MSE values, at various J/S ratios, deviates less, as 
compared to MFOM method. For example, from Table 6.6, we see that the average 
degradation computed, when MAML method is applied, extends from 48.68% up to 
78.5%, for -10 dB up to 10 dB J/S, while the MFOM method has a range from 54.36% 
up to 93.68%. 

The previous observations for the performance of the denoising methods, are also 
verified if we consider the improvement of the MSE values computed from the denoised 
data as compared to the ones obtained by the direct cross-correlation. 

In Table 6.8 the average percentage improvement of the MSE is summarized for 
the three data sets at each jamming-to-signal ratio. Note that the average SNR 
improvement values of the HP-ADS data set with the sampling interval 92.59 (nsec) 
appears to be inconsistent, if we consider the observations we have made so far. This is 
not true, because in this case, the computed values of the mean squared error are based 
upon a smaller range of SNRs (we are limited in SNR according to the J/S ratio). The 
improvement in this case must be quantified by taking the smaller SNR range into 
account. 
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Jamming- 

Average Percentage Improvement % 

to-signal 

Sampling Interval 

Sampling Interval 

Sampling Interval 

ratio 

370.370 (nsec) 

185.185 (nsec) 

92.59 (nsec) 

(dB) 


Modified 


Modified 


Modified 


MAML 

4 th Order 

MAML 

4 th Order 

MAML 

4 th Order 



Moment 


Moment 


Moment 






(20 dB down to 

(20 dB down to 






6dB SNR) 

6dB SNR) 

10 

84 

44.20 

82.40 

45.82 

89.8 

7.88 






(20 dB down to 

(20 dB down to 






6 dB SNR) 

6 dB SNR) 

5 

86.50 

58.71 

83.65 

47.05 

92 

37.27 






(20 dB down to 

(20 dB down to 






3 dB SNR) 

3 dB SNR) 

0 

88.60 

77.94 

86.10 

54.24 

85.34 

41.13 






(20 dB down to 

(20 dB down to 






0 dB SNR) 

0 dB SNR) 

-10 

90.60 

78.74 

86.90 

64.49 

77.46 

51.28 


Table 6.8. Average percentage improvement in MSE when Jamming is present 
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VII. SUMMARY, CONCLUSIONS AND FUTURE WORK 


A. SUMMARY AND CONCLUSIONS 

This thesis is focused on the improvements made possible by denoising, in 
conjunction with wavelet processing, for estimating the Time Difference Of Arrival 
(TDOA) of GSM signals received at two spatially separated sensors. 

Two wavelet based denoising methods are proposed to improve the TDOA 
estimation: the Modified Approximate Maximum Likelihood (MAML), and the Modified 
Fourth Order Moment (MFOM) method. 

The performance of these methods was evaluated under unequal SNR conditions 
as well as when jamming was present. The GSM signals used were generated with the 
help of the HP-ADS software. The Mean Squared Error (MSE) was used as criteria of 
goodness (low MSE denotes a small error and hence good localization). 

The study indicates that the application of the denoising method improves the 
TDOA estimation, as compared to when no denoising is applied to the received signals. 
In particular, the MAML method allows a 70.4 to 91.37 percent improvement in the 
TDOA estimation, while the respective improvement obtained when the MFOM method 
is applied, ranges from 62 to 79 percent. Increasing the sampling interval over which 
signal samples are taken from the signal decreases the MSE. 

In addition, we observed that the MAML method provided the best results, and 
outperformed the MFOM method, in all cases. The low SNR limits in each channel, 
above which denoising is beneficial were obtained. The MAML method appeared to be 
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more robust in presence of Gaussian noise antijamming, while the MFOM method shows 
an undesirable “sensitivity” at the lower SNR levels (i.e., below -3 dB). 

B. FUTURE WORK 

Follow on work should examine the possibility for a modification of the proposed 
algorithms, so that the time-varying characteristics of the signal can be accounted for. 

In addition, the performance of the proposed denoising methods should be 
investigated under the presence of a fading environment, and/or when other jamming 
schemes are applied (which may be more effective or severe), than the “erfed” FM noise 
we used. 

Finally, one could further quantify the results of our work by relating the MSE 
improvement to the proper solution of the hyperbolic equations as presented in Chapter 

m. 
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APPENDIX A. SIMMULATION RESULTS OF THE COMBINED 

DENOISING METHODS 

As mentioned in Section IV.C, combinations of the proposed denoising methods 
were tried, see Figure 4.3, to obtain their performance in a noisy environment. Figures 
A.l, A.2, and A.3, present the results of our simulations when we apply the combined 
methods (MAML/MFOM and MFOM/MAML), on the three HP-ADS GSM signals, 
embedded in noise. We use the same SNR in both channels, in a range from 20 dB to -6 
dB, and again the MSE stands as criteria of goodness (i.e. low MSE denotes a small error 
and hence good localization). 

Each figure, one for each of the HP-ADS data sets, depicts the MSE, obtained 
using the combined methods, versus SNR, and also includes the MSE values obtained 
when each denoising method (MAML or MFOM), is applied individually. The reader 
can easily compare the performance of the combined methods in this way. 

We have to point out that, from Figures A. 1(c), A.2(c), and A.3(c), one can see 
that the combination of the MAML/MFOM performs better than the MFOM alone. So 
we can use this combination instead of the MFOM. But this is not what we did. The 
reason is obvious if we take a look at Figure A.4 (a-d). This figure presents the results 
for MSE versus SNR, in the range of 20 dB to -6 dB for the HP-ADS data with a 
sampling interval T s =92.59 (nsec). It is clear that reasonable results using the 
MAML/MFOM combination are available for SNR values greater than -3 dB. It appears 
that the combined method can not perform well for low SNR values when the sampling 
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MSE MSE 


interval decreases, while the individually applied MFOM permits the extraction of 


reasonable MSE values, even in SNRs lower than -3 dB. 


(a) (b) 




(c) (d) 




Figure A.l. Combined denoising methods as compared to the MAML method and the 
MFOM method for the HP-ADS data with T s =370.37 (nsec) 
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MSE MSE 


(a) (b) 




(c) 



(d) 



SNR (dB) 


Figure A.2. Combined denoising methods as compared to the MAML method and the 

MFOM method for the HP-ADS data with T s =185.185 (nsec) 
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MSE MSE 


(a) 


(b) 




(c) (d) 




Figure A.3. Combined denoising methods as compared to the MAML method and the 

MFOM method for the HP-ADS data with T s =92.59 (nsec) 
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Figure A.4. Combined denoising methods as compared to the MAML method and the 
MFOM method for the HP-ADS data with T s =92.59 (nsec) at a range of SNR values 

from 20 dB to -6 dB 
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MSE MSE 


APPENDIX B. PERFORMANCE IN A JAMMING SITUATION 
(SCENARIO WITH 200 KHZ JAMMING BANDWIDTH) 

A. RESULTS OF JAMMING FOR HP-ADS DATA WITH T s =370.37 (NSEC) 


(a) J/S ratio 10 dB (b) J/Sratio5dB 




(c) J/S ratio 0 dB 




-5 0 5 10 15 20 

SNR 


Figure B.l. MSE versus SNR at different jamming-to-signal (J/S) ratios (data with 

T s =370.370 (nsec)) 
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Jamming-to-signal ratio 

(dB) 

Average Percentage Degradation % 

for HP-ADS data with sampling Interval 370.37(nsec) 

MAML method 

Modified 4 th Order Moment 

method 

10 

76.01 

96.86 

5 

73.12 

92.60 

0 

69.32 

81.10 

-10 

54.34 

56.60 


Table B. 1. Average percentage degradation in MSE at various jamming-to-signal 

ratios (J/S) 


As we can see from Figure B.l and Table B.l, the conclusions discussed in 
Section VI.B are also applicable to the current situation. In this case where a more severe 
jamming situation is present, the denoising methods provide improved results. The 
MAML method outperforms the MFOM method and appears to be more robust. The 
MFOM method, as is evident from Figure B.l (a-d), becomes more vulnerable as the J/S 
ratio increases. Of course the average degradation in MSE values, in this jamming 
scenario, is greater for both methods when compared to the corresponding values in 
Table 6.5, but this is expected since the jamming conditions are more intense. 
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B. RESULTS OF JAMMING FOR HP-ADS DATA WITH T s =185.185 (NSEC) 


(a) J/S ratio 5 dB (b) J/S ratio 0 dB 




(c) J/S ratio-10 dB 



Figure B.2. MSE versus SNR at different jamming-to-signal (J/S) ratios (Data with 

T s = 185.185 (nsec)) 


71 























Jamming-to-signal ratio 

(dB) 

Average Percentage Degradation % 

for HP-ADS data with sampling Interval 185.185 (nsec) 

MAML method 

Modified 4 th Order Moment 

method 

10 

N/A 

N/A 

5 

82.73 

98.94 

0 

82.62 

93.34 

-10 

67.60 

69.78 


Table B.2. Average percentage degradation in MSE at various jamming-to-signal 

ratios (J/S) 


C. RESULTS OF JAMMING FOR HP-ADS DATA WITH T s =92.59 (NSEC) 

In this jamming situation the application of the denoising methods did not yield 
reasonable results. The MSE values obtained were very high and the denoising was not 
beneficial. Only in the case that the jamming-to-signal (J/S) ratio is below -10 dB, we 
are able to obtain meaningful MSE values. The plot for -10 dB J/S ratio, can be viewed 
in the Figure B.3. 

Summarizing, we can say that the proposed denoising methods are beneficial even 
in a severe jamming environment. 
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Figure B.3. MSE versus SNR for HP-ADS data with sampling interval T s =92.59 

(nsec), at J/S ratio -10 dB 
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APPENDIX C. MATLAB CODES 

A. MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD 
(MAML) 


S-'k'k-k'k'k'k-k'k'k'k-k-k'k'k-k'k'k'k-k-k'k'k-k'k'k'k'k-k'k'k-k-k'k'k-k-k'k'k-k'k'k'k-k-k'k'k-k'k'k'k-k-k'k'k-k'k'k'k-k'k 

o 

%**********************************************************■*■■*■ 

o 

% DenoiseMAML: Modified Approximate Maximum Likelihood delay 
% estimation via orthogonal wavelet transform. 

% In this function we perform the wavelet 

% transform of the signal, compute the weights, 

% modify the decomposition coefficients, and 

% recombine the modified coeficients, using the 

% inverse wavelet transform, to obtain the 

% denoised signal. 


SYNTAX: y 


DenoiseMAML(xn,yn) 


INPUT: xn 
yn 


Received signal 
Received signal 


from first receiver 
from second receiver 


OUTPUT:y = Denoised signal 


% SUB_FUNC: None 
% Written by Spiros Mantis 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k 

^•k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


function y=DenoiseMAML(xn,yn) ; 

xyn=xcorr(xn,yn,'biased'); 
[sigmas b]=max(xyn); 
rx=xcorr(xn, 'biased') ; 
maxx=rx(length(xn)); 
ry=xcorr(yn,'biased'); 
maxy=ry(length(yn)); 
sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 

nx=wmaxlev(length(xn),'db4');; 

[cx lx]=wavedec(xn,nx,'db4'); 


dxc= [ ] ; 
for i=l:nx 

d=detcoef(cx,lx,i); 
dl=length(d) ; 
sigmad=(1/dl)*sum(d. A 2) ; 
sigmasd=sigmad-sigmanl; 


METHOD 
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sigmaa(i)=(2 A (nx-1))*sigmasd; 
if sigmasd<=0 
wd=0; 
else 

wd=sigmasd/(sigmanl*sigman2+sigmasd*(sigmanl+sigman2)); 

end 

dc=wd*d; 
dxc=[dc dxc]; 

end 

a=appcoef(cx,lx,'db4',nx); 
al=length(a); 
aux=sum(sigmaa) ; 
sigmasa=((2 A nx))*sigmas-aux; 
if sigmasa<=0 
wa=0; 
else 

wa=sigmasa/(sigmanl*sigman2+sigmasa*(sigmanl+sigman2)); 

end 

ac=wa*a; 
dxc=[ac dxc]; 

xd=waverec(dxc,lx,'db4'); 
y=xd; 
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S- icic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
o 

S^ici(*icic*i?icici(*icicici?icicic*ic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

% MAML_main: This is a test program for the Modified 
% Approximate Maximum Likelihood method. 

% In this program we used GSM signals 

% which were generated with the help of 

% the HP-ADS software. 

o, 

o 

% SYNTAX: MAML_main 

o. 

o 

% INPUT: None 

o, 

o 

% OUTPUT: Mean Square Error (MSE) versus SNR 

Q, 

O 

% SUB_FUN: DenoiseMAML.m 

Q, 

O 

% Written by Spiros Mantis 

o 

S^ic*ici(ic*ici(ici(ici(ici(ici(ic*ici(ic*ici('k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

clear all 

data_370_setl 
s=transpose(s2); 
sl=length(s); 

pow=(1/si)*sum(abs(s) . A 2) ; 


K=100 % number of realizations 

rand( 1 seed',40); 
f=150*rand(l,K) ; 
delay=floor(f); 

delay(1:K/2)=-delay(1:K/2); % delay is between 

% -150 to +150 samples 


n=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 

SNR=10. A (n./lO) ; 

for k=l:length(SNR) 

oran=SNR(k) 
for i=l:K 

x=[zeros (1,200) s zeros (1,224)]; 

y=[zeros(1,200+delay(i)) s zeros(1,224-delay(i))]; 
randn('state',2*(i+j)); 

noil_real=sqrt(pow/(2*oran))*randn(1, 1024) ; 
randn('state',3*(i + j) ) ; 

noil_imag=sqrt(pow/(2*oran))*randn(1, 1024) ; 
randn('state',4*(i + j)) ; 

noi2_real=sqrt(pow/(2*oran))*randn(1, 1024) ; 
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randn(’state',5*(i+j)); 

noi2_imag=sqrt(pow/(2*oran))*randn(1, 1024) ; 

noil=noil_real +j *noil_imag; 
noi2=noi2_real +j *noi2_imag; 

xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 


% TDOA calculation with xcorr(without denoising) 
xy=xcorr(xn,yn); 

[al bl]=max(xy); 

erl(i)=delay(i)-(bl-1024); 

% Denoising using the MAML method 
x_real=DenoiseMAML(real(xn) , real(yn)) ; 
y_real=DenoiseMAML(real(yn) , real(xn)) ; 
x_imag=DenoiseMAML(imag(xn),imag(yn)); 
y_imag=DenoiseMAML(imag(yn),imag(xn)); 

X=x_real +j *x_imag; 

Y=y_real +j *y_imag; 

% TDOA calculation with xcorr (with denoising) 
XY=xcorr(X,Y); 

[a2 b2]=max(XY); 

er8(i)=delay(i)-(b2-1024) ; 


end 

errorMAML(k) = (1/length(er8))*sum(er8. A 2) ; 
errorlaMAML(k) = (1/length(er1))*sum(er1. A 2) ; 

Hmaml(k,:)=er8; 

Hlamaml(k,:)=erl; 


end 


figure (1) 

k=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 

plot(k,errorlaMAML(1:14),'o',k,errorMAML(1:14),'x',k,errorlaMAML(1:14), 

k,errorMAML(1:14)) 

legend('Direct CC','MAML') 

title('Denoise MAML method for HP-ADS data 370.370 (nsec)') 

ylabel('MSE') 

xlabel('SNR') 

figure (2) 

plot(1:2047,xy) 

title ('Correlation Function of x and y signals with noise') 
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figure (3) 
plot(1:2047,XY) 

title('Correlation Function of x and y signals without noise') 


save errorMAML; 
save errorlaMAML; 

save Hmaml; 
save Hlamaml; 


79 



B. MODIFIED FOURTH ORDER MOMENT METHOD (MFOM) 


; DenoiseMFOM: Modified Fourth Order Moment delay estimation 
; via orthogonal wavelet transform. 

; In this function we perform the wavelet 

> transform of the signal, modify the 

; decomposition coefficients by thresholding, 

; and recombine the modified coeficients, 

; using the inverse wavelet transform, to 

; obtain the denoised signal. 

> SYNTAX: y = DenoiseMFOM(xn,yn) 

; INPUT: xn = Received signal from first receiver 
: yn = Received signal from second receiver 


OUTPUT: y = Denoised signal Xn based on Yn statistics 


% SUB_FUNC: None 
% Written by Spiros Mantis 

^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

^•k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


function y=DenoiseMFOM(xn,yn) ; 

xyn=xcorr(xn,yn,'biased'); 
[sigmas b]=max(xyn); 
rx=xcorr(xn, 'biased') ; 
maxx=rx(length(xn)); 
ry=xcorr(yn,'biased'); 
maxy=ry(length(yn)); 
sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 

lamdax=3.l*sigmanl A 2; 
nx=wmaxlev(length(xn),'db4'); 
[cx lx]=wavedec(xn,nx,'db4'); 

dxc= [ ] ; 
for i=l:nx 

d=detcoef(cx,lx,i); 
dl=length(d); 

for j=l:dl 

if d(j) A 2<1.76*sigmanl 
d(j)=0; 
else 

d ( j)=d(j); 

end 

end 
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A=(1/dl)*sum(d. A 4); 


if A<lamdax 

dc=zeros(1,dl); 
else 

dc=d; 

end 

dxc=[dc dxc]; 

end 

a=appcoef(cx,lx,'db4',nx); 
al=length(a); 
for m=l:al 

if a(m) A 2<1.76*sigmanl 
a (m) =0; 
else 

a (m) =a (m) ; 

end 

end 

A=(1/al)*sum(a. A 4) ; 

if A<lamdax 

ac=zeros(1,al); 
else 
ac=a; 
end 

dxc=[ac dxc]; 

xd=waverec(dxc,lx,'db4'); 

y=xd; 
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^-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

S^ici(*icici(*i(ici(*icic**icicic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

% MFOM_main: This is a test program for the Modified 
% Fourth Order Moment method. 

% In this program we used GSM signals 

% which were generated with the help of 

% the HP-ADS software. 

o. 

o 

o, 

o 

% SYNTAX: MFOM_main 

o. 

o 

% INPUT: None 

o, 

o 

% OUTPUT: Mean Square Error (MSE) versus SNR 

Q, 

O 

% SUB_FUNC: DenoiseMFOM.m 

Q, 

O 

% Written by Spiros Mantis 

S^ic*ici(icicici(ici(ici(ici(ici(i(*ici(ic*ic'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

S^ic*ici(ic*ici(ici(ici(ici(ici(ic*ici(ic*'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

clear all 

data_370_setl 
s=transpose(s2); 
sl=length(s); 

pow=(1/si)*sum(abs(s) . A 2) ; 


K=100 % number of realizations 

rand('seed',40); 
f=150*rand(l,K) ; 
delay=floor(f); 

delay(1:K/2)=-delay(1:K/2); % delay is between 

% -150 to +150 


n=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 

SNR=10. A (n./lO) ; 

for k=l:length(SNR) 

oran=SNR(k) 
for i=l:K 

x=[zeros (1,200) s zeros (1,224)]; 

y=[zeros(1,200+delay(i)) s zeros(1,224-delay(i))]; 
randn('state',2*(i+j)); 

noil_real=sqrt(pow/(2*oran))*randn(1, 1024) ; 
randn('state',3*(i + j) ) ; 

noil_imag=sqrt(pow/(2*oran))*randn(1, 1024) ; 
randn('state',4*(i + j)) ; 

noi2_real=sqrt(pow/(2*oran))*randn(1, 1024) ; 
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randn(’state',5*(i+j)); 

noi2_imag=sqrt(pow/(2*oran))*randn(1, 1024) ; 

noil=noil_real +j *noil_imag; 
noi2=noi2_real +j *noi2_imag; 

xn=x+noil; % x + noise 
yn=y+noi2; % y + noise 


% TDOA calculation with xcorr(without denoising) 
xy=xcorr(xn,yn); 

[al bl]=max(real(xy) ) ; 

erl(i)=delay(i)-(bl-1024) ; 

% Denoising using the MFOM method 
X_real=DenoiseMFOM(real(xn),real(yn)); 
X_imag=DenoiseMFOM(imag(xn),imag(yn)) ; 
Y_real=DenoiseMFOM(real(yn) , real(xn)) ; 
Y_imag=DenoiseMFOM(imag(yn),imag(xn)); 

X=X_real +j *X_imag; 

Y=Y_real +j *Y_imag; 

% TDOA calculation with xcorr (with denoising) 
XY=xcorr(X,Y); 

[a2 b2]=max(real(XY) ) ; 

erlO(i)=delay(i)-(b2-1024); 


end 

errorMFOM(k)=(1/length(er10))*sum(er10. A 2); 
errorlaMFOM(k)=(1/length(er1))*sum(er1. A 2); 

Hmf om(k, :)=er10; 

HMFOM(k,:)=erl; 


end 


figure(1) 

k=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 

plot(k,errorlaMFOM(1:14),'o',k,errorMFOM(1:14),'x',k,errorlaMFOM(1:14), 

k,errorMFOM(1:14)) 

legend('Direct CC','MFOM') 

title ('Denoise MFOM method for HP-ADS data 370.370 (nsec) ') 

ylabel('MSE') 

xlabel('SNR') 

figure (2) 

plot(1:2047,xy) 

title('Correlation Function of x and y signals with noise') 
figure (3) 
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plot(1:2047,XY) 

title ('Correlation Function of x and y 


signals without 


noise ' ) 


save errorMFOM; 
save errorlaMFOM; 

save Hmfom; 
save HMFOM; 
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C. JAMMING 


Q'iri^iri^iri^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o 

% Jammer: Genaration of the noise jamming waveform. 

% This programm generates the "Erfed" FM 

% Noise, used as jamming waveform in our 

% simmulations. 

o, 

o 

% SYNTAX: Jammer 

Q, 

O 

% INPUT: None 

o. 

o 

% OUTPUT: 1) "Erfed" FM noise jamming waveform 

% 2) Power spectrum of the "Erfed" FM noise 

% 3) Power spectrum of the HP-ADS signal 

Q, 

O 

% Written by: Spiros Mantis 

o, 

o 

o 

o 

clear all 

% Load the data generated from the HP-ADS 

data_370_setl; 

s=transpose(s2); 

sl=length(s); 

pow=(1/si)*sum(abs(s). A 2); 
rand('seed',40) 

% Define the desired SNR 
SNR=10 A (50/10); 
oran=SNR; 

y=[zeros (1,200) s zeros (1,224)]; 

% Generation of the noisy signal 
randn('state',2*(i+j)) 

noise_real=sqrt(pow/(2*oran))*randn(1, 1024) ; 
randn('state',3*(i + j) ) 

noise_imag=sqrt(pow/(2*oran))*randn(1, 1024) ; 
noise=noise_real +j *noise_imag; 
yn=y+noise; 

% FFT of the noisy signal 
S_spect=fft(yn) ; 

S_spectl=S_spect.*conj(S_spect); 

S=fftshift(S_spectl); 

Sl=S/max(S); 

figure (1) 

plot(10*logl0 (SI)) 

% Define the desired signal-to-jamming ratio 
oranj = 10 A (0/10); 
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% Define Time Vector -tmax to tmax 
N=1024; 

tmax=1024/(2*(2.7*1CT6)); 
t=2 *tmax*(-N/2:N/2-l)/N; 

% Sample Frequency fs 
fs=N/(2*tmax); 

% Generate Gaussian Noise Vector 
randn('state',6*(i+j)) 
n=randn(1,N); 
figure (2) 

subplot(2,2,1);plot(t,n);grid; 

xlabel('Time');ylabel('Noise Voltage'); 

title ('Gaussian Noise Voltage'); 

% Generate 'Erfed' Noise 
yfm=erf(n); 

% Shape Noise with Elliptical Filter 
[b,a]=e11ip(6,0.5,0.5,0.009) ; 
v=filter(b,a,n) ; 
verf=filter(b,a,yfm); 

% Define VCO with k=300 KHz/v 
f0=0;k=300000; 
f=f0+k*v; 
fl=f0+k*verf; 

x=sqrt(pow/(2*oranj))*cos(2 *pi*f.*t)+j*sqrt(pow/(2*oranj))*sin(2*pi*f.* 
t) ; 

xerf=sqrt(pow/(oranj))*cos(2*pi*fl.*t)+sqrt(pow/(oranj))*sin(2*pi*f1.*t 

); 


% Find Power Spectrum for v(t) 

V=fft(v); 

Pv=V.*conj(V) ; 

Pv=fftshift(Pv) ; 

% Find Power Spectrum for x(t) and xerf(t) 
X=fft(x); 

XE=fft(xerf); 

Px=X.*conj(X); 

Pxe=XE.*conj(XE); 

Px=fftshift(Px); 

Pxe=fftshift(Pxe); 


% Normalize Amplitude 

Pxl=Px/max(Px);Pvl=Pv/max(Pv);Pxel=Pxe/max(Pxe); 
Pxdbm=10*logl0(Pxl) ; 

Pvdbm=10*logl0(Pvl) ; 

Pxedbm=10*logl0(Pxel) ; 

% Select frequency axis 
fsam=fs*(-1024/2:1024/2-1)/1024; 
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% Plot control noise spectrum 
subplot(2,2,2); 
plot(fsam,Pvdbm);grid; 
xlabel('Hz');ylabel('DBM'); 
title('Noise Spectrum'); 

% Plot 'Erfed' FM noise spectrum 

subplot(2,2,3) ; 

plot(fsam,Pxedbm,'r');grid; 

axis([-1.2*10^6,1.2*10^6,-50,30]) ; 

xlabel('Hz') 

ylabel('DB'); 

title('Erfed FM Noise Spectrum'); 

% Plot Erfed FM noise waveform 

subplot(2,2,4); 

plot(1:1024,real(xerf) ) 

xlabel('Samples'); 

title('Jamming waveform');grid 

% Generation of the NOISE+Jamming+SIGNAL 
ynl=y+noise+xerf; 

% FFT of the SIGNAL+NOISE+JAMMING 
S_spectjam=fft(ynl) ; 

S_spectljam=S_spectjam.*conj(S_spectjam); 
Sjam=fftshift(S_spectljam); 

S1j am=S j am/max(S j am) ; 

figure (3) 

subplot (3,1,2) 

plot (1:1024,real(xerf)) 

title('Erfed waveform') 

subplot(3,1,1) 

plot (1:1024,real(yn)) 

title('Noisy signal') 

subplot(3,1,3) 

plot(1:1024,real(ynl)) 

title('Signal + Noise + Jamming') 
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^'k'k'k-k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k 


o 

o, 

o 

% Jamming_main: This is a test program the MAML method or the 
% MFOM method, when jamming is present. 

% In this program, we used the HP-ADS GSM signals. 

% We also used the "Erfed" FM noise jamming 

% waveform, which is added only in one channel. 

o, 

o 

% SYNTAX: Jamming_main 

o. 

o 

% INPUT: None 

o, 

o 

% OUTPUT: Mean Square Error (MSE) versus SNR 

Q, 

O 

% SUB_FUN: DenoiseMAML.m 
% DenoiseMFOM.m 

o. 

o 

% Written by Spiros Mantis 

%********************************************************■*■■*■■*■■*■■*■■*■■*■ 

o 

o 

clear all 

data_370_setl 
s=transpose(s2); 
sl=length(s); 

pow=(1/si)*sum(abs(s) . A 2) ; 

K=100 % number of realizations 

rand('seed',40); 
f=150*rand(l,K) ; 
delay=floor(f); 

delay(1:K/2)=-delay(1:K/2); % delay is between 

% ^150 to +150 


% Define the desired signal-to-jamming ratios 
Jammer=[-10 -5 0 10]; 

for m=l:4 

n=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 
SNR=10. A (n./lO); 
oranj=10 A (Jammer(m)/10); 

% Define Time Vector -tmax to tmax 
N=1024; 

tmax=N/(2*(10.8*10 A 6)); 
t=2*tmax*(-N/2:N/2-l)/N; 

% Sample Frequency fs 
fs=N/(2*tmax); 

% Generate Gaussian Noise Vector 
randn('state',6*(i+j)); 
noi=randn(1,1024); 

% Generate 'Erfed' Noise 
yfm=erf(noi); 

% Shape Noise with Elliptical Filter 

88 



[bfm,afm]=e11ip(6,0.5,0.5,0.04); 
v=filter(bfm,afm,noi); 
verf=filter(bfm,afm,yfm) ; 

% Define VCO with k=500 KHz/v 
f0=0;kfm=500000; 
fm=f0+kfm*v; 
fl=f0+kfm*verf; 

xfm=sqrt(pow/(2 *oranj))*cos(2 *pi*fm.*t)+j *sqrt(pow/(2*oranj))*sin(2*pi* 

fm.*t) ; 

xerf=sqrt(pow/(oranj))*cos(2*pi*fl.*t)+sqrt(pow/(oranj))*sin(2*pi*fl.*t 

); 


for k=l:length(SNR) 

oran=SNR(k); 
n (k) 

Jammer(m) 

for i=l:K 

x=[zeros(1,200) s zeros(1,224)]; 
y=[zeros(1,200+delay(i)) s zeros (1,224- 
delay (i)) ] ; 

randn('state',2*(i+j)); 
noil_real=sqrt(pow/(2*oran))*randn (1, 1024); 

randn('state',3*(i+j)); 
noil_imag=sqrt(pow/(2*oran))*randn(1,1024); 

randn('state',4* (i + j)); 
noi2_real=sqrt(pow/(2*oran))*randn (1, 1024); 

randn('state',5* (i + j)); 
noi2_imag=sqrt(pow/(2*oran))*randn(1,1024); 
noil=noil_real +j *noil_imag; 
noi2=noi2_real +j *noi2_imag; 
xn=x+noil+xerf; % x + noise 

yn=y+noi2; % y + noise 

% TDOA calculation with xcorr( without 
denoising) 
xy=xcorr(xn,yn) ; 

[al bl]=max(real(xy)) ; 
erl (i)=delay(i)-(bl-1024) ; 

% Denoise using MAML method 
x_real=DenoiseMAML(real(xn) , real(yn)) ; 
y_real=DenoiseMAML(real(yn),real(xn)) ; 
x_imag=DenoiseMAML(imag(xn) , imag(yn)); 
y_imag=DenoiseMAML(imag(yn),imag(xn)); 
X=x_real +j *x_imag; 

Y=y_real +j *y_imag; 

% TDOA calculation with xcorr(with 
denoising) 

XY=xcorr(X,Y) ; 

[a2 b2]=max(real(XY)); 
er8(i)=delay(i)-(b2-1024); 


end 
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error8JAM(k) = (1/length(er8))*sum(er8. A 2) ; 
errorlaJAM(k) = (1/length(er1))*sum(er1. A 2) ; 

H8JAM(k,:)=er8; 

HlaJAM(k,:)=erl; 

end 
figure(m) 

k=[20 18 17 16 15 14 13 12 9 6 3 0 -3 -6]; 

plot(k,errorlaJAM(1:14),'o',k,error8JAM(1:14),'x',k,errorlaJAM(1:14),k, 
error8JAM(1:14)) 

legend('Direct CC','MAML') 

title('MAML method with Jamming for HP-ADS data 370.370 
(nsec) ' ) 
ylabel('MSE') 
xlabel('SNR') 

ERROR8JAM(m,:)=error8JAM; 

ERRORlaJAM(m,:)=errorlaJAM; 


end 

figure (5) 

plot(k,ERROR8JAM(1,:),'r',k,ERROR8JAM(2,:),'b',k,ERROR8JAM(3,:),'g',k,E 
RROR8JAM(4,:),'c') 

legend('-lOdB S/J ratio','-5dB S/J ratio','0dB S/J ratio','10dB S/J 
ratio') 

title('MAML method with Jamming for HP-ADS data 370.370 (nsec)') 
xlabel('SNR') 
ylabel('MSE') 
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